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Abstract 

The amorphous aluminium silicate (Al203)2(Si02) [AS2] is investi- 
gated by means of large scale molecular dynamics computer simulations. 
We consider fully equilibrated melts in the temperature range 6100 K 
> T > 2300 K as well as glass configurations that were obtained from 
cooling runs from T — 2300 K to 300 K with a cooling rate of about 
10 12 K/s. Already at temperatures as high as 4000 K, most of the Al and 
Si atoms are four-fold coordinated by oxygen atoms. Thus, the structure 
of AS2 is that of a disordered tetrahedral network. The packing of AIO4 
tetrahedra is very different from that of S1O4 tetrahedra in that Al is in- 
volved with a relatively high probability in small-membered rings and in 
triclusters in which an O atom is surrounded by four cations. We find as 
typical configurations two-membered rings with two Al atoms in which 
the shared O atoms form a tricluster. On larger length scales, the system 
shows a microphase separation in which the Al-rich network structure 
percolates through the S1O2 network. The latter structure gives rise to a 
prepeak in the static structure factor at a wavenumber q = 0.5 A -1 . The 
comparison of experimental X-ray data with the results from the simu- 
lation shows a good agreement for the structure function. The diffusion 
dynamics in AS2 is found to be much faster than in SiC>2. We show that 
the self-diffusion constants for O and Al are very similar and that they 
are by a factor of 2-3 larger than the one for Si. 
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1 Introduction 



Amorphous mixtures of SiC>2 with other oxides such as Na20 or AI2O3 are of 
fundamental interest in geosciences and glass technology These systems 
exhibit many of the physical phenomena that one encounters in multicompo- 
nent melts: E.g. in sodium silicates the sodium ions are much more mobile 
than the silicon and oxygen atoms which leads to the property of ion conduc- 
tance at low temperature (Hj. Or in the mixture Si02~Al2 03 a miscibility gap 
emerges which is already present at temperatures slightly below 2000 K 
Some recent computer simulations have shown that in order to understand these 
properties of silicate melts one needs an accurate knowledge of their micro- 
scopic structure 1 \ For instance, pure amorphous silica forms a network of 
corner-shared SiC>4 tetrahedra and it is known from Molecular Dynamics (MD) 
computer simulations that diffusive motions are dominated by the existence of 
defects such as SiC>3 and SiOs units in the network (see, e.g., Refs. 8 9). An- 
other example is the sodium diffusion in sodium silicates: Molecular dynamics 
simulation studies have demonstrated that the sodium diffusion can be under- 
stood in terms of a motion through a channel network which is embedded in 
the SiC>2 matrix |10lllllH2lll3j , and that the characteristic length scale of this 
channel network is directly related to a prepeak in the static structure factor 
at q = 0.95 A -1 |14U15j the existence of which has been verified in a recent 
neutron scattering experiment Despite this progress in our understanding 
of the microscopic structure of these materials, many properties of silicates are 
still understood only poorly. E.g. the influence of a second network former, such 
as Al, on the structural and dynamical properties, or on the phase diagram are 
still not very clear. Therefore, the goal of the present paper is to shed some 
light on this matter. 

Aluminium silicate glasses have been investigated by means of different ex- 
perimental techniques such as Nuclear Magnetic Resonance (NMR) |17lll8ll*H?l 
E0IEI1E21E3 , IR and Raman spectroscopy, and X-ray scattering [2"ll2"5ll2l)] . 
One of the main issues in these experiments has been the analysis of the local 
structure around the aluminium atoms. This is of special importance for under- 
standing the chemical ordering in aluminium silicates since in principle the Al 3+ 
ions do need a different environment of O 2- ions than the Si + ions in order to 
provide local charge neutrality. Noteworthy are two peculiarities that have been 
found in the experimental studies that distinguish the local oxygen environment 
of an Al from that of a Si atom: In systems with a high AI2O3 content such as 
mullitc 3(Al 2 03)2(Si02) there seems to be a relatively large amount of five- and 
six-fold coordinated Al atoms in addition to AIO4 units. And secondly, NMR 
experiments found evidence for the existence of a high amount of so-called 
triclusters, i.e. structural units where an oxygen atom is surrounded by three 
cations (whereby at least one of them is an aluminium atom) fl9.20.21.2*21l23j. 
The possibility of such triclusters has been very recently confirmed in ab initio, 
molecular orbital calculations |27j . The presence of triclusters has been also dis- 
cussed in the context of extensive viscosity measurements of Na20-Al203~Si02 
liquids EEHU]. 
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As mentioned before, pure silica forms a network of corner-shared SiC>4 units 
and thus the O atoms are two-fold coordinated by Si atoms. The appearance 
of triclusters in amorphous aluminium silicates is an indication for the different 
local ordering of the different cationic species Al and Si. The difference in the 
local structure around the Al atoms from that around the Si atoms is accompa- 
nied by the tendency towards a metastable liquid-liquid phase separation below 
» 1900 K for an A1 2 3 content between about 10 to 50 mol% HEE]. Th e sys- 
tem we consider in the following is approximately in the centre of the demixing 
region: By means of MD simulations we investigate an aluminium silicate melt 
with 33 mol% AI2O3, i.e. (Al203)2(Si02), and in the following we will denote 
this system by AS2. 

It is one of the merits of MD simulations that one can study at the same time 
the structure on local and medium length scales in conjunction with dynamic 
properties. The duration of the MD runs is of course restricted to time scales 
in the ns range and so we are not able to approach the experimentally expected 
liquid-liquid coexistence line very closely. However, as we demonstrate in the 
following, the simulations do shed light onto precursors of phase separation in 
the microscopic structure. 

The rest of the paper is organized as follows: In the next section we give the 
details of the model and the simulation procedure. Then we present the results 
for the structure and diffusion in AS2, and finally we summarize and discuss 
them in the last section. 

2 Model and details of the simulation 

A potential that has often been used recently to study the properties of amor- 
phous silica is the so-called BKS potential 3QI that has been developed by van 
Beest, Kramer, and van Santen by means of ab initio calculations. Although 
it is a simple pair potential, it has been shown to reproduce very well many 
static and dynamic properties of amorphous silica [3"Tll3"2"ll3"3'l [34 8,35, 331EZIEH1 
I39II4Q|IH1EE21|43II44| . An extension of the BKS potential that allows to consider 
also mixtures of silica with other oxides such as Na20 and AI2O3 was proposed 
by Kramer et al. gSJ. As demonstrated recently |j"o1ll71lrll4iairolOll1fln^ . 
this potential gives also a quite realistic description of the static and dynamic 
properties of sodium silicates. 

In this work we use the latter potential to investigate the aluminium silicate 
melt (Al2C>3)2(Si02) [AS2]. The functional form of the potential is as follows: 

/ \ qaqpe 2 . , f r> \ Ca/3 
tPaf3[r) = — \- A a pexp(-B al3 r) — , (1) 

where a, j3 € {Si, Al, O}. Here r is the distance between an ion of species 
a and an ion of species (3. The values of the parameters {A a ^ 1 B a p 1 C a p\ 
that were calculated by ab initio methods are Asio = 18003.7572 eV, Aaio = 
8566.5434 eV, A 00 = 1388.7730 eV, B Si p = 4.87318 A" 1 , 5 A io = 4.66222 A -1 , 
Boo = 2.76 A"\ Csio = 133.5381 eVA 6 , Caio = 73.0913 eVA 6 , and Coo = 
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175.0 eVA 6 (for the Si-Si, Si-Al and Al-Al interactions the latter parameters 
are all set to zero) |3B|. 

From Eq. it becomes obvious that at small distances the potential be- 
tween the Al (or Si) and the O atoms goes to minus infinity (since the coefficients 
C a p are positive), i.e. it becomes unphysical. Therefore we have modified the 
potential at short distances by substituting it with a simple parabola: 

eV / ? \ ? q 

<Pmo = -19.9508 eV + 15.5 (r - 1.16576 A) r< 1.1658 A. (2) 

A 

We emphasize however, that this change hardly affects the properties of the 
system since at intermediate and low temperatures the probability of finding a 
distance below 1.17 A between an Al and an O particle is very small. Note that 
we use similar modifications as the one given by Eq. J2J) also for the Si-0 and 
0-0 interactions. Details can be found in Ref. |31j . 

In the long-ranged Coulomb-part the charges q a e (e: charge of an electron) 
are not the bare ionic charges of ions of type a but are considered to be effective 
charges. Unfortunately in Ref. |45j . the g Q 's were fixed such that systems like 
AS2 are not neutral: go = —1-2, qsi = 2.4, and qai = 1.9. Only with the addi- 
tional component phosphorus the system recovers charge neutrality according 
to the parameter sets in Ref. 05] • The same problem arises for sodium sili- 
cates if one uses the Kramer potential. In order to overcome this problem we 
follow here the same strategy that we have already applied in the case of the 
potential for the sodium silicates {46114148) : We set the charge of aluminium to 
qps\e = 1.8e such that a stoichiometric aluminium silicate system as AS2 is neu- 
tral. Then we add a short-ranged potential that compensates for this change 
at short distances. The final form of the potential is given by 

^ = ^ + ^^[l-(l-^Ai)(l-^Ai)]e(r c -r) , (3) 
with O being the Heavisidc function, with <jsi — 2.4, qo = —1.2, and 

(r c - rf 



QAl{r) = <?ai 1 + In 



C 



lA 2 + (r c - r) 2 



Here the parameters r c — 6.0 A, C = 0.0653, and d = 2.0 A 2 are chosen. These 
parameters guarantee that the potential given by Eq. J2J is indeed very similar 
to the original one as illustrated in Fig.Q] For r < 2.2 A the modified potential 
is almost the same as the original one. Finally we mention that the presence of 
the exponential term in Eq. Q makes that the charge (7ai is a smooth function 
of r. 

The equations of motion were integrated with the velocity form of the Ver- 
let algorithm using a time step of 1.6 fs. The Coulomb forces were calculated 
by the standard Ewald summation technique. At each temperature we equili- 
brated the system first in the NVT ensemble by coupling it to a stochastic heat 
bath. Thereby, the equilibration time exceeded the structural relaxation time 
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for the slowest component (i.e. silicon), i.e. dynamic density-density correla- 
tion functions for silicon at the wave-vector q = 1.7 A -1 (corresponding to the 
characteristic length scale of the tetrahedral network, see below) have decayed 
to zero in that time [5JJ . After equilibration the heat bath was switched off, and 
we started production runs in the microcanonical ensemble. During all runs the 
density was fixed to the experimental value at 300 K, p = 2.60 g/cm 3 [52|. The 
number of particles was 1408 (iVgj = 256, TVai = 256, No = 898) and hence the 
linear dimension of the cubic simulation box was L — 26.347 A. The tempera- 
tures that we have investigated are 6100 K, 4700 K, 4000 K, 3580 K, 3250 K, 
3000 K, 2750 K, 2600 K, 2480 K, 2380 K, and 2300 K. The production runs at 
the lowest temperature lasted over 6.9 ns real time corresponding to 4.2 million 
time steps. In addition we did cooling runs from 2300 K to K with a constant 
cooling rate 7 = 1.42 • 10 12 K/s. At each temperature we did five completely 
independent runs in order to improve the statistics. 

3 Results 

In this section we present the results of our simulations. First, we study in detail 
the local structure of the network in AS2 and show then the consequences of 
the local chemical ordering for length scales that go beyond distances of nearest 
and next-nearest atoms. Finally, we discuss the behavior of the self diffusion 
constants. 

Before we start the discussion of structural quantities we present in Fig. [5] 
the pressure p as a function of temperature. From this graph we recognize that 
in the temperature range considered the pressure varies between 1 and 5 GPa. 
This strong variation is related to the fact that we have made our simulation 
at constant volume. It is, however, reassuring that at ambient temperature the 
pressure is not exceedingly high, thus showing that the potential is quite reliable 
with respect to the pressure. Furthermore we point out that it is unlikely that 
for a system like AS 2 the structure changes significantly in the pressure range 
that we have here. Therefore we expect that the results presented here are very 
similar to the ones that one would get in a constant pressure simulation. 

It is remarkable that p(T) exhibits a local minimum around T = 2500 K 
which would correspond to a density maximum in a constant pressure simula- 
tion. So our model predicts for AS2 a density anomaly which is a well-known 
feature in amorphous silica Note that for T < 2300 K the system is no 

longer in equilibrium but in the glass phase. From the fact that in this temper- 
ature range p(T) shows basically a linear variation with temperature, we thus 
can conclude that anharmonic effects are not important. 

Quantities that are well suited to characterize the local structure of atomic 
systems are the partial pair correlation functions g a p (r) which are proportional 
to the probability of finding a particle of type (3 at a distance r from a particle of 
type a. The definition of g a p(r) can be found in standard textbooks The 
six different g a p(r) of our system are shown in Fig.|3|for the three temperatures 
T = 4000 K, 2300 K, and 300 K. The sharp first peak around r = 1.605 A in 
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3sio(V) reflects the strong covalent nature of the Si-0 bond. We also can note 
that in the first minimum of gsio(f), which is around 2.3 A, the function is 
basically zero even at 4000 K. Thus this feature allows for a natural definition 
of nearest neighbors of a silicon atom, and below we will make use of this fact. 
Moreover, there is a gap between 2 A and 3 A which is due to the chemical 
ordering in the network of AS2: In between a silicon atom and a second nearest 
oxygen neighbor there must be always another silicon (or aluminium) atom. The 
functions gAio(?~) look very similar to <?sio( r ) at the corresponding temperature, 
but the former are less pronounced in that the first peak is slightly broader 
and it has a smaller amplitude. Furthermore, the position of the first peak in 
<7aio(?") is now at the slightly higher value r — 1.66 A. This is in agreement 
with experiments and ab initio simulations of similar systems, although slightly 
higher values between 1.71 A and 1.77 A have been reported |55115!)]. 

Significant differences are found between 3aiai(V) and <?siSi(f) : At 300 K the 
first peak in gAiA\{r) splits up in two peaks at n = 2.59 A and r2 = 3.16 A 
whereas in gsiSi( r ) one finds only a single peak at about 3.12 A, i.e. a value 
very similar to r 2 . Also in <7siAi(?") one finds a shoulder around n. We will see 
below that the feature around r% is due to the presence of two-membered rings. 
Fig-Et shows a comparison of <?oo( r ) for AS2 with that for SiC>2 at T = 300 K 
(the latter was taken from a recent simulation study, for details see Ref. [S]). 
The main difference is that the first peak in the function for AS2 is shifted to 
larger distances which stems from the fact that the length of an Al-0 bond is 
slightly larger than the length of a Si-0 bond. 

In pure SiC>2 a disordered tetrahedral network is formed such that a silicon 
atom sits in the centre of each tetrahedron, whereby the oxygen atoms at the 
four corners of this tetrahedron are shared by the silicon atoms of the two 
neighboring tetrahedra (and thus each oxygen atom is two-fold coordinated by 
silicon atoms). It was shown in a MD simulation of a SiC>2 model |S]) that this 
local structure is essentially formed at temperatures as high as 3000 K since 
even at this temperature the percentage of defects (such as a Si atom that is 
five-fold coordinated by O atoms or an O atom that is three-fold coordinated 
by Si atoms) is smaller than 5 %. 

We demonstrate now that the structure of our AS2 model is far from a 
perfect tetrahedral network even at very low temperatures. To this end, we 
consider coordination number distributions P a /3(z) which give the probability 
that a particle of type a is surrounded by exactly z neighbors of type /3 within 
a distance r < r^f n (where rj^f n corresponds to the first minimum in g a p(r)). 
Fig. H shows P af} (z) for Si-O, Al-O, 0-(Si,Al), O-Si, and O-Al correlations 
(0-(Si,Al) means that one does not distinguish between Si and Al atoms) at 
the three temperatures 4000 K, 2750 K, and 300 K. The values used for rf^L are 
r^9 = 2.33 A, r£P = 2.54 A for T = 4000 K, r^P = 2.25 A, r^P = 2.42 A for 
T = 2750 K, and r^P = 2.20 A, r^-° = 2.40 A for T = 300 K. As we can infer 
from PsiO and Paio > at low temperatures most of the Si and Al atoms are four- 
fold coordinated by O atoms and thus also Al is integrated in the tetrahedral 
network structure. Note that also a recent MD simulation of a realistic model 
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of a pure AI2O3 melt finds a structure at low temperatures in which most of the 
Al atoms are four-fold coordinated by O atoms [ST] . 

That the packing of the AIO4 tetrahedra is nevertheless different from that 
of the SiC>4 tetrahedra is demonstrated in Figs. 0J: and^Ji. Only 70% of the O 
atoms are two-fold coordinated but around 30% of the O atoms are three-fold 
coordinated by (Si,Al) atoms, thus forming the so-called triclusters that have 
been identified in NMR experiments. It is remarkable that the percentage of 
such triclusters is nearly independent of temperature in our simulations. The 
distributions Posi and Poai show that the probability that an O atoms is three- 
fold coordinated by silicon atoms is very low (essentially zero at 300 K), whereas 
the probability for triclusters with three Al atoms is relatively high and increases 
even slightly with decreasing temperature. A closer inspection at T = 300 K 
shows the following cation composition of the triclusters: There are 1.4% with 
three Si, 12.4% with two Si and one Al, 47.3% with one Si and two Al, and 
38.9% with three Al. So most triclusters contain one Si and two Al atoms or 
three Al atoms. 

We have seen that at low temperatures most of the Al atoms are four-fold 
coordinated by oxygens. On the other hand there are two different characteristic 
length scales 7*1 and r 2 for the distance between nearest Al neighbors. These two 
length scales should be also reflected in the geometry of the AIO4 tetrahedra: 
Two connected tetrahedra for which the Al atoms are at a distance r\ from 
each other may have a different geometry from two connected tetrahedra where 
the two aluminium atoms in the centres are at distance around r 2 ■ Appropriate 
quantities to study the geometry of the tetrahedra are P a p^{9), the distribution 
for the O-Si-0 and O-Al-0 angles, which are shown in Fig.EJ In Posio(^) a 
single peak is observed which, by decreasing the temperature, becomes sharper 
and the location of the maximum moves to larger angles. At T — 300 K the 
maximum is at 9 = 108.2° which is close to the value for an ideal tetrahedron, 
9 = 109.47° (vertical lines). The behavior of -Poaio(^) is very different in that 
one finds at T = 300 K two peaks: 9 X = 85.8° and 9 2 = 109.8°. Such a 
bimodal distribution has been found also in computer simulations of free silica 
surfaces .38:; 39^ where one obtains two peaks in Posio(^) for the surface region 
at similar values for 9\ and 9 2 . And also in this case the bimodal distribution 
of angles is accompanied by two "nearest neighbor peaks" in gsisK 7 ")- The 
authors of Refs. [3*%ll3l?] have explained these features at free silica surfaces by 
the presence of two-membered rings, i.e. structural units where two tetrahedra 
share two oxygen atoms. The two tetrahedra that form the two-membered rings 
are deformed such that the distance between the two cations is at r\ and the 
O-cation-0 angles are at 9\ . (Here the two oxygens are of course the shared O 
atoms.) There is one major difference between the two-membered rings of free 
silica surfaces and the ones in AS2: In the former case they become very rare at 
low temperatures whereas in the case of AS2 their occurrence increases slightly 
with temperature. 

Apart from two-membered rings it is also straightforward to define rings of 
any size n in the network: One selects any cation (i.e. Si or Al) and two of 
its O neighbors. There are several possible paths how one can move through 
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the network structure from one cation-0 pair to the next one such that one 
starts from the pair with the first O atom and ends at a cation-0 pair with the 
second O atom. The shortest possible of such paths are called rings and the 
length of a ring is the number of cation-0 pair that it contains. In Fig. [S]the 
ring size distribution P(n) is shown for three different temperatures. Whereas 
P(n) changes significantly from 4000 K to 2300 K, only minor changes are 
observed from 2300 K to 300 K. However, we recognize essentially the same 
qualitative behavior at the three temperatures: most probable are rings with a 
size of n = 5 and there is a relatively large contribution of rings with n = 2 and 
n = 3. Note that the location of the maximum at n — 5 is different from that 
which one expects for pure silica at low temperatures where MD simulations 
have found the maximum at n — 6 (corresponding to the length of the rings in 
/3-cristobalite) [3TU58] . 

We have also determined the composition of small-membered rings at T = 
300 K yielding for the rings with n = 2 that 72.8% consist of two Al atoms, 
25.0% contain one Al and one Si atom, and there are 2.2% with two Si atoms. 
For the three-membered rings we obtain the following numbers: 36.5% with 
three Al atoms, 45.3% with two Al and one Si atom, 15.5% with one Al and 
two Si atoms, and 2.7% with three Si atoms. 

One may ask whether the relatively high probability of two-membered rings 
in our AS model is related to the occurrence of the aforementioned triclusters. 
And indeed we have extracted from our data that at T = 300 K 96% percent of 
the oxygen atoms that are involved in a tricluster are also a member of a ring 
with n = 2. And this holds also the other way round: more than 95% of the O 
atoms that are members in any two-membered ring arc three-fold coordinated 
by cations and form thus a tricluster. So we find that the appearance of triclus- 
ters in our AS2 model is accompanied by the presence of two-membered rings 
which to our knowledge is a connection that has not been considered yet in the 
analysis of experiments. A typical local configuration with a two-membered 
Al-0 ring and two triclusters is illustrated by a schematic picture in Fig. [7| 

The results discussed so far have shown that the aluminium atoms do form 
tetrahedral units with oxygens but that the packing of these AIO4 units is very 
different from the one of the Si04 units. We now want to study whether the 
different local order around Al and Si atoms leads to structural features that are 
present on larger length scales. To this end, we consider partial static structure 
factors which are the Fourier transforms of the corresponding pair correlation 
functions. They are correlation functions of the number densities 



N, 





(5) 
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(depending on wave-vector q) and can be defined as follows |54) : 

Sap(q) = <Pa(q)p/3(-q)) 

= ^EEN ,( !'^- r ')^ ( 6 ) 

k=l 1=1 

with TV being the total number of particles, q is the absolute value of the wave- 
vector q, and f a p is equal to 0.5 if a =^ /3 and to 1.0 if a = 0. Fig. |H] shows 
S a p(q) at the three temperatures T = 4000 K, 2300 K, and 300 K. Similar to the 
results found for the radial distribution functions, the temperature dependence 
of S a p(q) is relatively weak in that the different peaks that are found at low 
temperatures are already present at a temperature as high as 4000 K. For q > 
2.3 A -1 the partial structure factors reflect length scales of nearest neighbors 
and their location corresponds approximately to the period of oscillations in the 
9a/3(r)- But we are now more interested in the features in S a p(q) at small q. 
The peaks at 1.6-1.7 A -1 are due to the order that arises from the tetrahedral 
network structure, i.e. from repeated AIO4 and SiC>4 units. And indeed the 
length 27r/(1.7 A ) = 3.7 A corresponds approximately to the spatial extent 
of two connected tetrahedra. Note that a prepeak around q = 1.7 A -1 in the 
static structure factor is also found in pure silica and in many other materials 
that form tetrahedral networks (e.g. see Ref. J3H and references therein). 

But from Fig. [S] we recognize that there is also an additional prepeak at 
q = 0.5 A -1 in the S a p(q) for the Si-Si, Al-Al, and Si-Al correlations. One 
possible explanation for this peak is the following one: As we have seen in our 
analysis of the local structure before, the local packing of AIO4 tetrahedra is 
significantly different from that of the SiC>4 tetrahedra. This may lead to a 
structure where an AIO4 tetrahedron prefers to be surrounded by other AIO4 
tetrahedra and thus Al rich regions are formed. So a structure is created where 
the Al atoms are not homogeneously distributed on the relatively large length 
scale I » 27r/(0.5 A~ ) w 12.6 A. If one considers only the Al atoms for instance, 
voids are formed with the spatial extent given by I in the regions where the Si 
atoms sit. The same holds of course true if one considers only the Si atoms. 
In the S a p(q) in which the O atoms are involved, almost no peak at 0.5 A -1 
is seen because the oxygen atoms are essentially homogeneously distributed in 
the system on the length scale I since they are both nearest neighbors of Si and 
Al atoms (the difference in the chemical ordering of the O atoms around the Al 
atoms from that around the Si atoms is only weakly pronounced on the length 
scale I). 

The structure of our AS2 model at T = 300 K is illustrated by a snapshot 
of the simulation box, Fig. (Note that the size of the shown atoms does not 
correspond to their actual size.) One can clearly see that Al rich regions are 
formed that percolate through the Si04 network. From this figure it is also 
visible that the packing of the AIO4 tetrahedra is denser due to the formation 
of more compact structural units such as small-membered rings. 
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If our interpretation of the prcpeak at 0.5 A -1 is correct then it should be 
present in a pronounced way in static concentration fluctuations. For a mix- 
ture with the two components A and B the static concentration-concentration 
structure factor S cc (q) can be easily computed from the partial structure factors 

S cc {q) = x|5 A a(?) + x 2 A S B B(q) - 2x A x B SAB(q) (7) 

with xa = Na/N and xb = Nb/N being the concentration of A and B particles, 
respectively. For our three-component system we may also use Eq. (0 in the 
following way: We do not distinguish between two of the three species such that 
for instance the A species combines Si and O particles and the B species is given 
by the Al particles. For this example the functions Saa(<z)j <Sab(<z)j an d Sbb(q) 
are then given by 

Saa(<?) = SsiSi + Soo 2<SsiO , 
Sab{q) = SsiAl + Saio , 

SbbG?) = Saiai j 

and from these functions a structure factor S cc (q) can be calculated by Eq. Q 
whereby Na = Nsi+No and N-q — Nai- The latter concentration-concentration 
structure factor is nothing else than the static autocorrelation function of the 
concentration density cai of the Al atoms |59llfi0| . In an analogous way one can 
calculate two other structure factors for AS2 that correspond to the autocorre- 
lation functions of the concentration densities for Si and O, respectively. Thus, 
in terms of the concentration densities 

Ca(q) = P«(q) ~ x a (psi(q) + PAi(q) + Po(<l)) a e {Si, Al, 0} (8) 

with x a = N a /N, the latter three concentration-concentration structure factors 
are defined by 

S CQCQ (9) = ^<c Q (q)c Q (-q)) . (9) 

More details on the S CaCa (q) can be found in Refs. |59U60) . 

The three functions S CaCa (q) are shown in Fig. ^] for T — 300 K. In all 
three functions a pronounced peak is seen around q = 2.72 A -1 corresponding 
to the length scale of nearest Si-Si, Al-Al etc. neighbors. This is reasonable 
since there should be strong concentration fluctuations on this length scale due 
to the chemical ordering in the tetrahedral network of AS2. Furthermore we 
do also find a well-pronounced peak around 0.5 A -1 in S CAlCAl (q) and S CSiCSi (q) 
which confirms our interpretation of this prepeak. Also in accordance with our 
interpretation is that no prepeak around 0.5 A -1 is seen in S coco (q) because 
the O atoms are ordered in a similar way around the cations and so in S coco (q) 
no distinction is made between Si rich and Al rich regions. 

Of course the question arises whether it is possible to observe the prepeak 
at 0.5 A -1 also in experiments, i.e. whether or not real AS2 shows a structure 
on the length scale of 12.6 A. Unfortunately in experiments such as neutron 
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scattering one does not have access to the partial structure factors for systems 
like AS2 (due to the lack of appropriate isotopes) and one can only measure 
a linear combination of the partial structure factors (in the case of neutron 
scattering the S a p(q) are weighted by the neutron scattering lengths). As it is 
shown elsewhere |51| . in such a quantity one can hardly identify the prepeak 
at 0.5 A -1 because at small wave- vectors the dominant contribution comes 
from Soo(q), a partial structure factor which does not show the prepeak. So 
it remains a challenge to the experimentalists to develop techniques with which 
one can verify the presence of the latter prepeak in AS2. 

However, in order to see how far our microscopic model is able to reproduce 
the structure of real AS2, we compare now the "reduced" total static X-ray 
scattering factor of AS2, q(Sx(q) — 1) ; as obtained from our simulation to an 
experimental result. Sx{q) can be calculated from the S a p(q) by weighting them 
with X-ray form factors: 

with a, (3 € {Si,Al, O}. The formfactors f a (s) depend on the wave-vector q 
via s — q/4:n. We have taken the f a (s) from Ref. Fig. ITT1 shows Sx{q) as 
calculated from our simulation by Eq. (|10|l in comparison to the experimental 
result of Morikawa et al. [23] for an aluminium silicate melt with 37.1 mol% 
AI2O3, i.e. a composition which is slightly different from that of AS2 (AS2 
contains 33.3 mol% AI2O3). As we recognize from Fig.llllthe agreement between 
simulation and experiment is particularly good for q < 2.3 A^ 1 and also a fair 
agreement is obtained for higher q. 

Finally, we want to turn our attention to the diffusion dynamics of the AS2 
melt. The self-diffusion constant D a for a particle of type a e {Si, Al, 0} can 
be calculated from the mean squared displacements (r^(t)) via the Einstein 
relation: 

A>=lim^M (11) 

In Fig. E| the three different D a for AS2 are plotted on a semi-logarithmic 
scale as a function of inverse temperature. Also included are Dsi and -Do of 
pure silica from a recent MD simulation [Hj. At T = 6100 K the D a {T) in AS2 
are very similar to those in Si02- However, upon decreasing the temperature, 
the dynamics in AS2 does not slow down as rapidly as the one of SiC>2. At T = 
2750 K the diffusion of all components in AS2 is about two orders of magnitude 
faster than in SiC>2. The self-diffusion constants in SiC>2 show a crossover from 
a power law behavior as predicted by the mode coupling theory (MCT) of the 
glass transition |S2] at high temperatures (see dashed lines in Fig. 1121) to an 
Arrhenius behavior at low temperatures whereby we have found the critical 
mode coupling temperature at 3330 K (for more details see Refs. jSlESI) • A 
similar analysis as the one in SiC>2 for AS2 is the subject of future work. As 
we see in Fig. ^] the self-diffusion in AS2 provides the very interesting case 
of a tetrahedral network where the dynamics of one of the cations, aluminium, 
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is slightly faster than that of oxygen whereas the diffusion of the other cation, 
silicon, is a factor of 2 or 3 slower than that of oxygen. 

In order to quantify further the difference in the temperature dependence of 
the different diffusion constants, we show in Fig.^|the ratios D$i/Do, D^i/Dxh 
Dai/ Dq for AS2 as well as Dgi/Do for SiC>2 as a function of temperature. 
Dsi/Do for SiC>2 varies only from 0.65 to 0.8 in the temperature range 3400 K < 
T < 6100 K, whereas it decreases rapidly below about 3400 K. The latter 
behavior can be explained by a change in the transport mechanism around the 
critical MCT temperature T c = 3330 K of our SiC>2 model (for more details see 
Refs. As mentioned before, in the case of AS2 it is an open question 

whether one can understand the diffusion dynamics by means of MCT. However, 
Dsi/Do and Dsi/Dai exhibit a very similar behavior as Dgi/Do in SiC-2. It is 
remarkable that Dai/ Do seems to approach a constant value for temperatures 
below about 3000 K, an issue that has to be clarified also in future studies. 

4 Summary and conclusions 

In this paper we have presented the results of large scale molecular dynamics 
computer simulations to investigate the structure and diffusion dynamics of 
the amorphous aluminium silicate melt (Al203)2(Si02) [AS2]. The microscopic 
interactions were described by a simple pair potential that has been determined 
by Kramer et al. |45j using ab initio calculations. 

It is one of the open questions in the literature on aluminium silicate melts 
and glasses how the local structure around an aluminium atom evolves in order 
to yield charge neutrality on a local scale. We find in the simulations of our 
AS2 model that, similar to the Si atoms, already at high temperatures most 
of the Al atoms are four-fold coordinated by O atoms. However, the packing 
of the AIO4 tetrahedra is denser than that of the SiC>4 tetrahcdra which is 
manifested especially in the presence of 3(Al,Si) triclusters. Evidence for such 
structural units in systems like AS2 has been given in NMR experiments |19I2()I 
l2Tll2*2ll2*3"| . Our simulation allows a detailed description of the geometry of the 
triclusters: We find that in most of them the three-fold coordinated O atoms 
are members of a two-membered ring (whereby these rings contain most likely 
two aluminium atoms). Furthermore, we find also a relatively high occurrence 
of three-membered rings in AS2 in which most likely at least two Al atoms 
participate. In contrast to the small mcmbcrcd rings in Si02, Al rings in AS2 
of length two or three become slightly more frequent if one cools down the 
system to low temperatures. 

Similar structural features as in our AS2 model have been found recently also 
in a MD simulation of pure amorphous aluminium oxide (AI2O3) • Also in 
that case a high number of small-membered rings was found, in particular also 
a relatively high number of two-membered rings. The geometry of the latter 
rings is slightly different from that in AS2 since in AI2O3 (at least as far as it 
is predicted in the simulation) it is not the presence of two edge-sharing AIO4 
tetrahedra that is very likely but either two edge-sharing AIO5 polyhedra or an 



12 



AIO4 tetrahedron sharing an edge with an AIO5 unit. However, the presence 
of two-membered rings in a model of AI2O3 and in a model of AS2 (based on 
completely different ab initio calculations) may indicate that two-membered Al 
rings are indeed a typical structural unit in melts and glasses containing AI2O3 
and thus, this finding is worth to be checked in experiments. 

The difference in the local environment of the Al and the Si atoms leads to 
the formation of Al rich regions such that a network of AIO4 units percolates 
through the SiC>4 structure. The characteristic length scale of the latter struc- 
ture is reflected by a prepeak in the static structure factor at q « 0.5 A -1 . This 
prepeak can be seen as the manifestation of a microphase separation and indeed, 
in accordance with this interpretation, it is a pronounced feature in static struc- 
ture factors of concentration densities as defined by Eq. (J7J) (see Fig. . The 
presence of a microphase separation in the melt structure is another prediction 
of our simulation that could be checked in experiments. Of course, as we have 
shown above, this is not a simple task. 

It is interesting that prepeaks with a physical origin as the one in our AS2 
model have been found also for other binary silicate melts, e.g. in sodium silicates 
at q = 0.95 A" 1 both in simulation and neutron scattering or in a calcium silicate 
glass at q = 1.3 A -1 by means of neutron scattering |63| . 

The critical point of the demixing transition in the experimental phase dia- 
gram of Al 2 3 -Si0 2 is at sa 1920 K for a mixture with w 30 mol% A1 2 3 [§]. 
Thus, if the model used in this work exhibits a similar phase diagram as real 
systems, AS2 is close to the critical composition and the lowest temperature 
used in this work above the glass transition, T = 2300 K, is only a few hundred 
Kelvin above the critical temperature. In this sense, the observed microphase 
separation can be seen as the precursor of a critical unmixing transition. In 
this context it is interesting to study how the dynamics is affected by a possi- 
ble interplay between structural relaxation and a critical slowing down. It has 
recently been shown that mode coupling theory (MCT) describes the slowing 
down due to structural relaxation in systems like Si0 2 or (Na 2 0)2(Si0 2 ) very 
well 8,35,48 (this holds above and around the critical temperature of MCT). 
In AS2 the case might occur where the scenario for structural relaxation as 
predicted by MCT interferes with the critical dynamics near a second order 
demixing transition. 

In order to shed light on these guesses it is necessary to determine the phase 
diagram of the model and this requires the use of techniques that are different 
from Molecular Dynamics: Well-suited in this case are Monte-Carlo simula- 
tions in the semi-grandcanonical ensemble in conjunction with multicanonical 
sampling These methods have also the advantage that they yield config- 
urations exactly on the coexistence line of the demixing transition which can 
be used as starting configurations for Molecular Dynamics simulations to study 
the microscopic structure and dynamics. The latter procedure has recently been 
successfully used for a symmetrical Lennard- Jones mixture |65| . Therefore it 
can be hoped that the application of this technique will also be useful for more 
complex systems like AS2 and thus will help to increase our understanding of 
this glass-former. 



13 



Acknowledgments: This work was supported by SCHOTT Glas. Generous 
grants of computing time on CRAY-T3E at the NIC Jiilich and the HLRS 
Stuttgart are gratefully acknowledged. 



References 

[1] See, e.g., J. F. Stcbbins, P. F. McMillan, and D. B. Dingwell (Eds.), Rev. in 
Mineralogy 32 (1995). 

[2] See, e.g., H. Bach and D. Krause (Eds.), Analysis of the Composition and 
Structure of Glass and Glass Ceramics (Springer, Berlin, 1999). 

[3] An overview on the recent research on ion-conducting mixtures can be 
found in K. L. Ngai, Proc. of the 4th International Discussion meeting in 
Complex Systems, J. Non-Cryst. Solids 307-310 (2002). 

[4] J. Hlavac, The Technology of Glass and Ceramics - An Introduction (El- 
sevier, Amsterdam, 1983). 

[5] H. Rawson, Properties and Applications of Glass (Elsevier, Amsterdam, 
1980). 

[6] J. F. MacDowell and G. H. Bcall, J. Am. Ceramic Soc. 52, 17 (1969). 

[7] A review of the simulation of silicate melts can be found in P. H. Poole, P. 
F. McMillan, and G. H. Wolf, Rev. in Mineralogy 32, 563 (1995). 

[8] J. Horbach and W. Kob, Phys. Rev. B 60, 3169 (1999). 

[9] C. Oligschleger, Phys. Rev. B 60, 3182 (1999). 

[10] M. D. Ingram, Phys. Chcm. Glasses 28, 215 (1987); Phil. Mag. B 60, 729 
(1989). 

[11] W. Smith, G. N. Greaves, and M. J. Gillan, J. Chem. Phys. 103, 3091 
(1995). 

[12] J. Oviedo and J. F. Sanz, Phys. Rev. B 58, 9047 (1998). 

[13] P. Jund, W. Kob, and R. Jullicn, Phys. Rev. B 64, 134303 (2001). 

[14] J. Horbach, W. Kob, and K. Binder, Phys. Rev. Lett. 88, 125502 (2002). 

[15] J. Horbach, W. Kob, and K. Binder, J. Phys.: Condons. Matter 15, S903 
(2003). 

[16] A. Meyer, H. Schobcr, and D. B. Dingwell, Europhys. Lett. 59, 708 (2002). 

[17] S. H. Risbud, R. J. Kirkpatrick, A. P. Taglialavore, and B. Montez, J. Am. 
Ceram. Soc. 70, C10 (1987). 



14 



[18] R. K. Sato, P. F. McMillan, P. Dennison, and R. Duprcc, J. Phys. Chem. 
95, 4483 (1991). 

[19] B. T. Poc, P. F. McMillan, B. Cote, D. Massiot, and J.-P. Coututes, J. 
Phys. Chem. 96, 8220 (1992); B. T. Poc, P. F. McMillan, C. A. Angell, 
and R. K. Sato, Chem. Geol. 96, 333 (1992). 

[20] R. H. Meinhold, R. C. T. Slade, and T. W. Dawies, Appl. Magn. Reson. 4, 
141 (1993). 

[21] M. Schmiicker and H. Schneider, Ber. Bunsenges. Phys. Chem. 100, 1550 
(1996). 

[22] M. Schmiicker, K. J. D. MacKenzie, H. Schneider, and R. Meinhold, J. 
Non-Cryst. Solids 217, 99 (1997). 

[23] M. Schmiicker, H. Schneider, K. J. D. MacKenzie, and M. Okuno, J. Eur. 
Ceramic Soc. 19, 99 (1999). 

[24] H. Morikawa, S.-I. Miwa, M. Miyake, F. Marum, and T. Sata, J. Am. 
Ceramic Soc. 65, 78 (1982). 

[25] P. F. McMillan and B. Piriou, J. Non-Cryst. Solids 53, 279 (1982). 

[26] M. Okuno, N. Zotov, M. Schmiicker, and H. Schneider, submitted to J. 
Non-Cryst. Solids (2003). 

[27] J. D. Kubicki and M. J. Toplis, Am. Mineral. 87, 668 (2002). 

[28] M. J. Toplis, D. B. Dingwell, K. U. Hess, and T. Lend, Am. Mineral. 82, 
979 (1997). 

[29] M. J. Toplis, D. B. Dingwell, and T. Lenci, Gcochim. Cosmochim. Acta 
61, 2605 (1997). 

[30] B. W. H. van Beest, G. J. Kramer, and R. A. van Santcn, Phys. Rev. Lett. 
64, 1955 (1990). 

[31] K. Vollmayr, W. Kob, and K. Binder, Phys. Rev. B 54, 15808 (1996). 

[32] J.-L. Barrat, J. Badro, and P. Gillct, Mol. Sim. 20, 17 (1997). 

[33] M. Hemmati and C. A. Angell, in Physics meets Geology, H. Aoki and R. 
Hcmley, Eds. (Cambridge University Press, Cambridge, 1998). 

[34] J. Horbach, W. Kob, and K. Binder, J. Phys. Chem. B 103, 4104 (1999). 

[35] J. Horbach and W. Kob, Phys. Rev. E 64, 041503 (2001). 

[36] J. Horbach, W. Kob, and K. Binder, Eur. Phys. J. B 19, 531 (2001). 

[37] P. Schcidler, W. Kob, A. Latz, J. Horbach, and K. Binder, Phys. Rev. B 
63, 104204 (2001). 



15 



[38] A. Roder, W. Kob, and K. Binder, J. Chcm. Phys. 114, 7602 (2001). 

[39] C. Mischler, W. Kob, and K. Binder, Comp. Phys. Comm. 147, 222 (2002); 
C. Mischler, Ph.D. Thesis, Universitat Mainz, 2002. 

[40] P. Jund and R. Jullicn, Phys. Rev. B 59, 13707 (1999). 

[41] M. Bcnoit, S. Ispas, P. Jund, and R. Jullicn, Eur. Phys. J. B 13, 631 (2000). 

[42] D. J. Lacks, Phys. Rev. Lett. 84, 4629 (2000). 

[43] I. Saika-Voivod, F. Sciortino, and P. H. Poole, Phys. Rev. E 63, 011202 
(2001). 

[44] M. S. Shell, P. G. Debenedetti, and A. Z. Panagiotopoulos, Phys. Rev. E 
66, 011202 (2002). 

[45] G. J. Kramer, A. J. M. de Man, and R. A. van Santen, J. Am. Chem. Soc. 
113, 6435 (1991). 

[46] J. Horbach and W. Kob, Phil. Mag. B 79, 1981 (1999). 

[47] J. Horbach, W. Kob, and K. Binder, Chcm. Gcol. 174, 174 (2001). 

[48] J. Horbach and W. Kob, J. Phys.: Condens. Matter 14, 9237 (2002). 

[49] S. Ispas, M. Bcnoit, P. Jund, and R. Jullicn, Phys. Rev. B 64, 214206 
(2001). 

[50] A. Winkler, J. Horbach, W. Kob, and K. Binder, to be published. 
[51] A. Winkler, Ph.D. Thesis, Universitat Mainz, 2002. 

[52] O. V. Mazurin, M. V. Streltsina, and T. P. Shvaiko-Shvaikovskaya, Hand- 
book of Glass Data, Part A: Silica Glass and Binary Silicate Glasses (El- 
sevier, Amsterdam, 1983). 

[53] R. Briickncr, J. Non-Cryst. Solids 5, 123 (1970). 

[54] J. -P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic, 
London, 1986). 

[55] M. Benoit, S. Ispas, and M. E. Tuckcrman, Phys. Rev. B 64, 224205 (2001). 

[56] R. D. Shannon, Acta Crystallogr. Sect. A 32, 751 (1976). 

[57] G. Gutierrez and B. Johansson, Phys. Rev. B 65, 104202 (2002). 

[58] J. P. Rino, I. Ebbsjo, R. K. Kalia, A. Nakano, and P. Vashishta, Phys. Rev. 
B 47, 3053 (1993). 

[59] J. Bletry, Z. Naturforsch. 31a, 960 (1976). 



16 



[60] J. Horbach, submitted to J. Non-Cryst. Solids (2003). 

[61] International Union of Crystallography (Ed.), International Tables for X- 
ray Crystallography (Kynoch Press, 1974), vol. 4, chap. 2.2, pp. 71-99. 

[62] W. Gotze and L. Sjogren, Rep. Prog. Phys. 55, 241 (1992); W. Gotze, J. 
Phys.: Condens. Matter 10, Al (1999). 

[63] P. H. Gaskcll, M. C. Eckersley, A. C. Barnes, and P. Chicux, Nature 350, 
675 (1991). 

[64] D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in 
Statistical Physics (Cambridge University Press, Cambridge, 2000). 

[65] S. K. Das, J. Horbach, and K. Binder, J. Chem. Phys. 119 (2003). 



17 



30.0 



20.0 



„ 10.0 



> 



4>. 



Al-Si 



4>. 



AI-AI 



Kramer et al. 
mod. potential 




0.0 1.0 2.0 3.0 4.0 

r[A] 



5.0 6.0 



Figure 1: Al-Si, Al-Al and Al-0 interaction potentials as proposed by Kramer 
et al. |45j (dashed lines) and in the modified version (solid lines) according to 
Eqs. 
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Figure 2: The pressure as a function of temperature. Note that the points for 
T < 2300 K are sampled during runs where the system was cooled down from 
2300 K to K with a cooling rate of 7 = 1.42 • 10 12 K/s. 
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Figure 3: Partial pair correlation functions g a p{v) for the temperatures T = 
4000 K, 2300 K and 300 K. a) 5siS i(r), b) g S io(r), c) <?oo(r), d) jaiaiW, e) 
5aioM, f) ffSiAl(r)- The curve in c) for SiC>2 at T = 300 K is taken from Ref. (H|- 
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Figure 4: Distribution of coordination numbers P a p(z) for the temperatures 
T = 4000 K, 2300 K and 300 K, a) P Si o{z), b) P M o(z), c) P (Si,Ai)(«), d) 
PoSi{z) and Poai(z)- The three curves that at z = 2 have the largest values 
correspond to Posi( 2 )- 
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Figure 5: Distributions of angles 9 for the temperatures T = 4000 K, 2300 K 
and 300 K, a) -Posio(^), b) -Poaio(#)- The vertical lines correspond to the ideal 
tetrahedron angle of 109.47°. 
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Figure 6: Distribution of rings P{n) for the temperatures T = 4000 K, 3000 K, 
2300 K and 300 K. See text for the definition of a ring. 
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Figure 7: Schematic picture of a typical local configuration in AS2. O atoms are 
the big grey spheres, Al atoms are the small grey spheres, and the light small 
sphere is a Si atom. The black lines between the atoms symbolize covalent Al-0 
and Si-0 bonds. The picture shows a two-membered Al-0 ring whereby the O 
atoms in this ring form triclusters. 



24 




Figure 8: Partial static structure factors S a p(q) for the temperatures T — 
4000 K, 2300 K and 300 K. a) S SiSi (q), b) S Si0 (q), c) S 00 (q), d) S A1A1 (q), 
c) S A \o{q), f) S siA i(q)- 
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Figure 9: Snapshot of AS2 at T = 300 K. The big white spheres are the sili- 
con atoms, the big black spheres are the aluminium atoms and the small black 
spheres are the oxygen atoms. Note that the size of the spheres does not corre- 
spond to the actual size of the atoms. 
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Figure 11: "Reduced" X-ray scattering factor q(Sx(q) — 1) as calculated from 
the simulation using Eq. (|10fl (solid line) in comparison to the experimental 
result by Morikawa et al. [21] (dashed line). 
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Figure 12: Arrhenius- plot of the diffusion constants of silicon, aluminium and 
oxygen: Comparison of the systems AS2 and SiC>2 (the data for SiC>2 is taken 
from Ref. The dashed lines are power law fits to the SiC>2 data (see text 
for details). 
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Figure 13: Temperature dependence of the indicated ratios of the diffusion 
constants D a /D (j (a, (3 G {Si, Al, O}). 



30 



